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Abstract. We present the results of a double analysis of the ionizing cluster in NGC 588, a giant HII region (GHR) 
in the outskirts of the nearby galaxy M33. For this purpose, we obtained ground based long-slit spectroscopy and 
combined it with archival ground based and space borne imaging and spectroscopy, in the wavelength range 1100- 
9800 A. A first modeling of the cluster was performed using integrated properties, such as the spectral energy 
distribution (SED), broad band colors, nebular emission H/3 equivalent width, the main ultraviolet resonance lines, 
and the presence of Wolf-Rayet star features. By applying standard assumptions about the initial mass function 
(IMF), we were unable to fit satisfactorily these observational data. This contradictory result led us to carry out 
a second modeling, based on a resolved photometric analysis of individual stars in Hubble Space Telescope (HST) 
images, by means of finding the best fit isochrone in color-magnitude diagrams (CMD), and assigning a theoretical 
SED to each individual star. The overall SED of the cluster, obtained by integrating the individual stellar SEDs, 
is found to fit better the observed SED than the best solution found through the integrated first analysis, but 
at a significantly later stage of evolution of the cluster of 4.2 Myr, as obtained from the best fit to the CMD. A 
comparative analysis of both methods traces the different results to the effects of statistical fluctuations in the 
up per end of the I MF, w hich are significant in NGC 588, with a computed cluster mass of 5600 Mq, as predicted 
bv lCervifio et all i2002h . We discuss the results in terms of the strong influence of the few most massive stars, 
six in the case of NGC 588, that dominate the overall SED and, in particular, the ionizing far ultraviolet range 
beyond the Lyman limit. 

Key words. Stars: evolution - Stars: luminosity function, mass function - (Stars:) Hertzsprung-Russel (HR) and 
C-M diagrams - Stars: Wolf-Rayet - ISM: individual objects: NGC 588 



■ 1. Introduction Inversely, modeling accurately a GHR from its emission 

. . . . spectrum requires a good knowledge of the ionizing spec- 

Despite their small number in comparison with Sun- like , c -, i • i • , , , , , 

1 .... trum of its source, which is one of the elements that control 

stars, massive stars are most influential in their host galax- , . , . . .. , , r ,, • • 

' . b the temperature and ionization state of the gas. 1ms lon- 

ies. They are responsible for the release in the interstellar •• , ■ » ii A i u j ■ _j 

J y , lzmg spectrum is unobservable, and can only be derived 

medium of most of the mechanical energy, metals such as r j i * j. n 

6 - r ' from a model of the stellar source. 

oxygen, neon, sulfur, etc. and possibly carbon, and hard 

radiation. However, their formation is still poorly under- The stellar content of a GHR can be determined in two 
stood. It is generally admitted that massive stars form wa y s: either h y observing the individual stars, or mak- 
mainly in instantaneous starbursts, with a power -law ini- m § assumptions regarding the distribution of their ini- 
tial mass functionjlMFlwith universal slope (e.g.. lMassevl tial masses - The flrst solutlon can be used to stu(b / thc 
2003), although!^ GSi estimates that this slope may cluster itself via colo r-magnitude diagrams (CMDs, e.g. 
comprise a natural scatter of a few tenths between differ- i Malamuth et al. || 1996[) or to synthesize an i onizing spec- 
ent clusters. Most massive stars are observed and studied t0 be input t0 P hotoionization models fielano et alj 
in very massive clusters embedded in giant HII regions I2QQ2D- However, such a procedure requires quite high reso- 
(GHRs), in which the ionized gas provides information lutlon and sensitivity for extragalactic objects, which can 
on the intensity of the (unobservable) stellar Lyman con- 111 S eneral be reached onl y with the HST ' In the vast ma " 
tinuum, via the equivalent width of the nebular H/3 line. of other cases ' onl y inte g rated Properties of the clus- 

ter can be used, and an IMF, general ly a continuous power 



ter can be used, and an inlf , generally a continuous power 
law, must be supposed. Nevertheless. lCervino et all l(2002l l 



Send offprint requests to: L. Jamet, e-m< 
Luc . Jamet@obspm . f r showed that if the total initial mass of a cluster is less than 



2 



L. Jamet et al.: The ionizing cluster of NGC 588 



~ 10 4 Mq, statistical fluctuations around the mean IMF 
can strongly affect the main diagnostics of the cluster, 
like the H/3 equivalent width W(R 0), as well as the de - 
termined ionizing spectrum (see also lCerviho et al1l2003|h 
The effects of these fluctuations need to be tested with 
observational data. 

In this article, we show that the cluster embedded 
in NGC 588 nicely illustrates the issue of IMF statisti- 
cal sampling. Section [3 describes all the data sets used, 
spectroscopic and photometric, ground based and space 
borne, either obtained by us or retrieved from different 
archives. Further data processing is detailed in Section 03 
Section ^describes a first modeling of the star cluster, us- 
ing a "classical" approach, aimed at fitting its integrated 
spectral properties with a model assuming an analytical 
IMF. In Section [5] we introduce the importance of mod- 
eling the effects of statistical fluctuations. In Section 
we measure the fluxes of the individual stars detected in 
HST images, and use the photometric results to analyze 
the stellar content of the cluster, by assigning a model SED 
to each individual star. The results from the two methods 
of modeling are significantly discrepant, and in Section 
we explore how this discrepancy can be attributed to the 
role of the few most massive stars, whose number is heav- 
ily affected by statistical fluctuations in the case of the 
low mass cluster ionizing NGC 588. 



2. Observational data and standard reduction 

We have performed ground based optical spectroscopic 
observations of NGC 588 and have retrieved the avail- 
able archival spectroscopic and imaging data sets from 
the Hubble Space Telescope (HST), the International 
Ultraviolet Explorer (IUE), and the Isaac Newton Group 
(ING). Table ^ contains a summary of the different data 
sets used in this study. 

2.1. Optical spectra 

Optical spectra were obtained in 1998 August 27/28 with 
the 3.5m telescope at the Centro Astronomico Hispano 
Aleman (CAHA, Calar Alto, Almeria). The observations 
were made with the dual beam TWIN spectrograph, us- 
ing the blue and red arms simultaneously, with the beam- 
splitter at 5500 A, and SITe-CCDs of 2000 x 800 15/x pix- 
els. Four exposures of 1800 s each were taken with the 
blue arm in the range 3595-5225 A, while two 1800 s ex- 
posures in the range 5505-7695 A and two 1800 s in the 
range 7605-9805 A were taken in the red arm (hereafter, 
the B, Rl and R2 ranges, respectively). The gratings T7 
in the blue and T4 in the red give dispersions of 0.81 
and 1.08 A/pixel in second and first spectral order, re- 
spectively. A slit, oriented at PA= 121° (Fig. QJ, of size 
240" x (1.2 ± 0.2)" provided a spectral resolution, simi- 
lar in both arms, of FWHM~ 2 A, as measured in the 
background sky emission lines. 



Table 1. Journal of Observations. 
Spectroscopy 

CAHA (PA=121°, width=1.2", 1998-08-27/28) 
ranee : A (A) A A (A/pixel) exposure (s) 
B : 3595-5225 0.81 4 x 1800 



Rl: 5505-7695 


1.08 


2 x 1800 


R2: 7605-9805 


1.08 


2 x 1800 


HST/STIS (PA= 


-154°, width= 


=0.2", 2000-07-13) 


range (A) 


exposure (s) 


data set 


1447-3305 


1440 


O5CO05010 


1075-1775 


480 


O5CO05020 


1075-1775 


2820 


O5CO05030 


IUE (PA=102°, large aperture, 


1980-10-03) 


range 


exposure (s) 


data ID 


1150-1980 


18000 


SWP10273 


1850-3350 


18000 


LWR08943 


Imaging 


HST/WFPC2 (Proposal: 5384, 1994-10-29) 


Filter 


exposure (s) 


data set 


F170W 


2 x 350 


U2C60701T/702T 


F336W 


2 x 160 


U2C60703T/704T 


F439W 


2 x 180 


U2C60801T/802T 


F547M 


2 x 100 


U2C60803T/804T 


JKT (1990-08-20/21) 


Band 


Filter 


exposure (s) 


H/3 


4861/54 


2 x 1800 


Ha 


6563/53 


1800 


Ha continuum 


6834/51 


1300 



These spectra were reduced, following the stan- 
dard procedure, by means of the IRAF 1 package 
noao . twodspec . longslit. Photometric calibration was 
achieved by means of three standard stars, BD+28 4211, 
G191-B2B and GD71, observed during the same night and 
setup, but through a wide slit. We estimated the result- 
ing photometric accuracy to be 3% in color, except for 
the A < 4000 A range (residuals are ~ 5% intrinsic to 
this range, ~ 15% with respect to the global spectrum), 
where the sensitivity curve drops, and for the A > 8000 
A range (accuracy ~ 30%), which is strongly affected by 
telluric absorption bands. The absolute photometric level 
was determined with ~ 10% accuracy. The reduction in- 
cluded the subtraction of the background measured in the 
slit aside the nebula, and consequently, the light of the 
underlying stellar population of M33 was removed from 
our spectra. 

2.2. IUE spectra 

We retrieved two spectra fr om the final archive IUE Newly 
Extracted Spectra (INES 2 :lRodriguez-Pascual et al.ll999t 
ICassatella et al.ll2000t iGonzalez-Riestra et al.ll2000|h The 

1 IRAF is distributed by the National Optical Astronomy 
Observatory, which is operated by the Association of 
Universities for Research in Astronomy, Inc., under cooper- 
ative agreement with the National Science Foundation. 

2 The INES System has been developed by the ESA IUE 
Project at VILSPA. Data and access software are distributed 
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SWP and LWR spectra cover a total range of 1150-3350 
A, and were obtained in the low-resolution mode (~5 A) 
through the large aperture (9.5" x 22"; cf. Fig.QJ. 

The two IUE spectra were retrieved from the archive in 
two formats: the integrated, fully calibrated spectra, and 
the two-dimensional frames, linearized in wavelength and 
position. The slit position and mean angular resolution 
were known with an accuracy of a few arcseconds, a scale 
comparable to that of the densest spatial structures of 
the cluster, and we re-determined them more accurately 
following the procedure detailed in § 13.21 Furthermore, due 
to significant uncertainties in the photometry of the LWR 
spectrum, we limited our use of IUE data to the SWP 
spectrum. 

2.3. Ground-based images 

We downloaded a set of images of NGC 588 from the Isaac 
Newton Group of telescopes archive. These images were 
taken with the lm Jacobus Kapteyn Telescope (JKT) at 
the Observatorio del Roque de los Muchachos (La Palma) , 
through the filters Ha+[NII], H/3, and a continuum nar- 
row band filter close to Ha. 

Line emission only images were obtained by subtract- 
ing a version of the continuum image scaled with stars 
common to the three filter images. Then, we divided the 
Ha+[NII] image by the H/3 one, thus establishing an ex- 
tinction map of the nebula (the intensity of the [Nil] 
A6548+6583 doublet being small compared to Ha, accord- 
ing to our spectra). We observed no significant gradient of 
extinction across the object. We also registered the images 
with respect to the HST images, in order to locate the slits 
of the different spectrographs over the nebula (Fig.|2||. 



2.4. HST images 

Wc retrieved five pairs of WFPC2 images from the HST 
archive, taken through the filters F547M, F439W, F336W, 
F170W and F469N. The images were calibrated with 
the on-the-fiy calibration system. Each pair of exposures 
within a filter was combined with crrej in order to per- 
form cosmic ray rejection. 

The images we re corrected for filter contam ination, 
as recommended bv lMcMaster fc Whitmorel l)2002|) . These 
corrections consisted of increasing the measured fluxes by 
0.044 magnitude for filter F170W, 0.006, for F336W and 
0.002, for F439 W. We also accounted for charge-transfer 
efficiency loss (|Whitmore et al.| Il999|) . This was neces- 
sary because of the low background levels in our images 
(^0.3 to 0.5 DN). In the case of individual stellar pho- 
tometry (Section 16. Ill, we applied directly the formulae of 
I Whit more et alJ l)l999h to correct the number of counts 
in each star. In the case of integrated fluxes (Section 
we used an average correction for each filter based on the 
rough locus of the cluster on the images (X ~ 400 and 



mm 



WNLW 



WN * 



-5 

Aa (arcsec) 



Fig. 1. HST image of NGC 588 in filter F439W and the 
spectral slits. The peculiar stars (WNL, WN) are dis- 
cussed in Section [6J 



and maintained by INTA through the INES Principal Center 
at LAEFF. http://ines.laeff.esa.es/ 



Y ~ 450) and on the approximate counts found in the 
stars that dominate the total cluster flux. These correc- 
tions are of 0.037, 0.035, 0.036 and 0.044 magnitude for 
filters F547M, F439W, F336W and F170W, respectively. 
We neglected the effects of geometric distortions, because 
the cluster is close to the center of the images, and the pho- 
tometric effects of these distortions are small compared to 
other errors. 

We computed the barycentric wavelengths of the fil- 
ters. For this, we multiplied the spectral response of each 
filter by a few SEDs characteristic of young OB associa- 
tions, and calculated the corresponding barycentric wave- 
lengths. For F547M, F439W, F336W, F170W and F469N, 
we respectively found 5470, 4300, 3330, 1740 and 4690 A. 
These wavelengths were then used to compute the values 
of the extinction laws to be applied. 



2.5. STIS spectra 

The HST archive contains three STIS spectra through 
NGC 588, covering the two wavelength ranges 1447-3305 
A (with the grating G230L) and 1075-1775 A (with the 
grating G140L), the two spectra covering the latter range 
having a resolution of ~1.5 A. We retrieved these spectra 
with on-thc-fly calibration. The spectra pass only through 
the brightest star of NGC 588 (Fig.^l, and we used them 
to identify its spectral class C Section [6. 21 below) . 
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Fig. 2. Ha+[NII] image of NGC 588 and CAHA and IUE 
slits. The contours represent the Ha continuum. 

3. Additional spectroscopic and photometric 
processing 

3.1. Nebular extinction 

Once the spectra were reduced, we performed a self- 
consistent procedure to deredden the nebular lines, based 
on the Balmer lines Ha through H5. First, the emission 
line fluxes were corrected for foregrou nd Galactic extinc- 
tion, which amounts to Eb-v — 0.045 l|Burstein fc Heilesl 



19841). t he Ga lactic law being here the ensemble of the 



Seatonl lll979h ultraviolet law and of the optical law of 
NandvetlaiTl|l975j) (as shown in Seaton's paper). Then 
a reddening law was fitted to the thus corrected Balmer 
decrement, taking into account the photometric and mea- 
surement errors of these emission lines 3 . This process was 
performed for each spatial increment along the slit inde- 
pendently, thus obtaining the reddening along the slit. The 
residuals from the actual line ratios to the fitted law were 
then plotted along the slit, giving statistical fluctuations 
around a mean offset value (due to a small residual photo- 
metric calibration) except at those locations where the un- 
derlying stellar population has a significant contribution 
of Balmer lines in absorption. The equivalent width of the 
lines at these locations was thus computed and the cor- 
rection applied to the intensity of the emission lines, then 
computing the correct ed value of redd ening. Assuming the 
LMC extinction law l)Howarthlll983), a priori adequat e 
due to the metallicity of the object i Vilchez et al.lll98'jj) . 
we derived an Eb-v curve along the slit, and found it to 



be almost constant, with a value of 0.11 ± 0.02 around 
the stellar cluster. This result is virtually independent of 
the selected extinction law, since different laws predict 
nearly the same relative extinctions within the H(5-Ha 
wavelength range. 

3.2. Spatial registration between spectra and images 

We have registered precisely the position of the CAHA 
slit with respect to the HST images. This is achieved by 
convolving the HST images with the atmospheric point 
spread function (PSF) at the time of the ground based 
observations, and comparing the spatial flux distribu- 
tions. In order to compute the atmospheric PSF we mod- 
eled the spatial flux distribution of a point source in the 
2D spectral frames in the foll owing way. W e assumed 
the atmospheric PSF to be a Moffat (1969) function 4 , 
h(r) = (1 + (r/R c ) 2 )~P . We selected a standard star ob- 
served just before NGC 588, and examined its profile along 
the slit around 4400 A. This profile was clearly asymmet- 
rical, indicating image distortion inside the spectrograph. 
We found that the observed star profile is very well repro- 
duced by a model of the form H(x) = J^j=i a%M(x — c.;), 
where x is the spatial increment variable, Cj the spatial 
location of the peak of each component and its weight, 
and M(x) is the Moffat PSF integrated across the slit 
width. Thus, we assumed that the signal of the observed 
object was convolved by a Moffat function determined by 
turbulence, truncated by the slit, and split into three com- 
ponents within the spectrograph, before hitting the detec- 
tor. The same analysis near 6600 A and 8700 A led to 
similar models. Once the PSF model was established, we 
proceeded as follows. 

The blue 2D spectral frame was multiplied by the re- 
sponse of the HST filter F439W and the result integrated 
along the wavelength axis; in this way we obtained from 
the spectral frames the spatial flux distribution of NGC 
588 through the HST filter F439W. Then, we applied our 
image degradation model to the HST image for different 
values of the Moffat parameters R c and 3 (that tend to 
vary with time), placed a synthetic slit on the blurred im- 
age, and extracted the profile integrated across the slit. 
R c , 8 and the central position of the slit were fit by maxi- 
mizing the correlation between the profiles extracted from 
the spectrum and from the image. We found R c ~ 1.1" 
and 3 k, 1.5, corresponding to a seeing of FWHM~ 1.6 
arcsec. The PSF and the slit position were used to calcu- 
late the aperture throughputs used in Section 0] and El 

We also used the two-dimensional frames and the 
F170W HST/WFPC2 image to determine accurately the 
slit position and angular resolution of the IUE data, as- 
suming a Gaussian PSF and using the same technique as 
for the optical spectra. The resulting correlation is very 
good for a FWHM of 6.4", with however a strong pho- 
tometric mismatch between the data: the IUE flux inte- 



A value of t he electron temperature of 11000 K 
llVflchez et al.llT988f) was assumed for the theoretical Balmer 
ratio. 



4 This has been shown to be the b est analytical re presenta- 
tion of atmospheric turbulence (e.g., lTerebizr]|l99ll) 
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grated within the F170W band was found to be higher 
than the HST one by a factor 1.41. This mismatch be- 
tween IUE and H ST absolute calibrations has been found 
also for NGC 604 llOonzalez Dekado fc Perezl2000h . HST 
absolute photometry is accurate to a few %, so we decided 
not to correct the HST image, but to divide the IUE spec- 
tra by the factor 1.41. 

3.3. Nebular continuum subtraction 

To analyze the spectral properties of the cluster itself, we 
removed the nebular continu um from our spectr a. For this, 
we calculated it according to lOsterbrockl ljl974|) . including 
the HI free-free, free-bound and two-photon components 
as well as the Hel free-free and free- bound emissions, as - 
suming a gas temperature of 11000 K i Vilchez et alll988|k 
we reddened the computed nebular continuum to match 
the observations, normalized it with the nebular H/3 emis- 
sion line intensity in the spatial zone of the optical slit 
selected to study the cluster (knowing the H/3 equivalent 
width with respect to the nebular continuum), and re- 
moved it. 



3.4. Aperture throughput correction 

Once the optical and UV slit positions and angular PSFs 
were determined (see Section l3~2f . we calculated the aper- 
ture throughputs of the slits for the whole cluster, i.e. the 
fraction of light transmitted by these slits. For this, we 
computed the sky background of each image by averaging, 
at a given position, the pixels contained in a rectangular 
box and close enough to the median value of this box. The 
resulting map of the sky and of the nebula was removed 
from the original in order to keep only the stellar signal. 
The remaining stellar images were then convolved by the 
optical or UV PSFs, and in each band, we compared the 
total flux of the cluster to the one entering the slit. In 
optical, the aperture throughputs were found to be inde- 
pendent of the band (F170W, F336W, F439W or F547M), 
and we retained a mean value of 0.169 for a slit width 
of 1.2", valid in the whole optical range. This through- 
put was proportional to the chosen slit width, known to 
«20% accuracy. The uncertainty on the PSF led to an un- 
certainty of 10% on the throughput. For the IUE spectra, 
we found a throughput of 0.83 in the F170W filter. This 
means that we had to divide the CAHA spectra by 0.169, 
and the IUE spectra, by 0.83, in order to obtain the inte- 
grated spectrum of the whole cluster. We also calculated 
the aperture throughputs for nebular emission: 0.062 for 
the optical spectra, and 0.44 for IUE data. 

4. Integrated canonical modeling of the cluster 

The cluster was first modeled in the classical way, i.e. 
as a stellar population whose IMF is rigorously a con- 



tinuous power law dN/dm oc 



TOlo 



< m < 



Assuming instantaneous star formation, we used 



Starburst99 



il.l ll999n . with miow =lMp 



0' 



making use of standard mass-loss evolutionary t racks 
at metalli cities Z=0.004 llCharbonnel et alJ Il993|) and 
Z=0.008 llSchaerer et"aTl Il993h clos e to the metallic- 
ity of the gas l)Vflchez et al.l Il988h (the solar value 
is he re Z=0.02). We sel ec ted the atmosphere mod- 
els of iLeieune et all lll997h . iHillier fc Miller! lll998h and 
IPauldrach et alJ (l200lt). t he lat ter two ones having been 
compiled bv lSmith et alJ l)2002f) . The age r of the cluster, 
its metallicity Z, and the IMF slope a and upper limit m up 
were constrained with observables available for the cluster 
as a whole. These are: i) the UV and optical colors, inte- 
grated over the entire cluster; ii) the H/3 equivalent width 
TV (H/3) measured over the whole object (nebula and clus- 
ter) ; iii) the profiles of the UV stellar lines present in the 
IUE spectrum; iv) the presence or absence of the main 
Wolf-Rayet (WR) emission lines in the optical SED. 



4.1. SED colors and derived extinction 

The SED colors were used to determine the amount and 
the law of extinction of the stellar flux. Indeed, the color 
excess Eb-v derived from the Balmer lines of a nebula 
has been found not to apply systematical ly to its ioniz- 
ing source (e.g., iMas-Hesse & Kunthll999l . Furthermore, 
whereas different extinction laws generally do not dif- 
fer significantly from each other in the optical range 
and cannot be discriminated by the use of the nebular 
Balmer lines, they strongly vary in the UV counterpart. 
Fortunately, the intrinsic UV and optical colors of an OB 
association little depend on the age and IMF of the latter, 
and can be used to characterize the extinction, provided 
that one knows approximately the main cha racteristics of 
this association. IMas-Hesse fc Kunthl l|l999f) found an age 
of 2.8 Myr for a very flat IMF (a — 1) with an upper 
mass cutoff of 120 M , for the cluster of NGC 588. We 
computed SEDs for a range of ages, IMF slopes and up- 
per mass cutoffs around these values, and calculated the 
corresponding fluxes in the four bands F547M, F439W, 
F336W and F170W. Then, the colors of the cluster were 
measured on the HST images cleared of their diffuse back- 
ground. Table shows the measured colors and some 
model outputs. The differences of color between the obser- 
vations and the models were attributed to interstel lar red- 
dening. We tested three extincti on laws: Galact ic l|SeatorJ 
rmTfllNandv et al]|l975h. FMC ftWrthll 9sl and SMC 
ijPrevot et alJll984UBouchet et alJll985h 7 whose values at 
the barycentric wavelengths of the bands are summarized 
in Table [21 For any set of model colors, the best fit was 
obtained with the SMC law and E B -v = 0.08 ± 0.04 (the 
error bar includes the uncertainties in the observational 
data and the spread in the theoretical colors) . This value is 
compatible with the amount of extinction derived from the 
nebular Balmer lines, meaning that there is no indication 
of differential extinction between the stellar and the neb- 
ular fluxes. Since the nebular value, Eb-v = 0.11 ± 0.02, 
was determined more accurately, we have also adopted 
it for the cluster in the following. Considering these re- 
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Filter 


F547M 


F439W 


F336W 


F170W 


A (A) 


5470 


4300 


3330 


1740 


Obs. M547 - M x 


0.00 


0.84 


1.48 


3.03 




±0.03 


±0.03 


±0.04 


±0.05 


M547 - M x (a) 


0.00 


0.96 


1.85 


4.10 


M547 - M A (b) 


0.00 


0.87 


1.49 


3.22 


M547 - M A (c) 


0.00 


0.82 


1.19 


2.94 


/(A) (Galactic) 


3.21 


4.18 


5.06 


7.80 


/(A) (LMC) 


3.10 


4.13 


5.37 


9.25 


/(A) (SMC) 


2.72 


3.65 


4.65 


10.51 



Table 2. Observed and model colors with filter F547M as 
the zero point. The observed colors are here corrected for 
foreground Galactic extinction. The three models (a), (b) 
and (c) shown here stand for m up — 120 Mq, Z=0.004, 
and respectively, a = 1, 1 and 2.35 and r —2, 3 and 4 Myr. 
We also show the three tested extinction laws in the form 
/(A) = A(X)/E B -v- The color differences A(M y - M\) 
between the observations and a given model were to be 
well fit by a function of the form a + Eb-vJ(^) for the 
tested extinction law to be accepted. 



suits and the 800 kpc distance of M33 l|Lee et alJl2002|) . 
we derived the cluster luminosity in the F439W band: 
4.71xl0 35 ergs -1 A" 1 . 

4.2. Hf3 equivalent width 

This quantity is often used in modeling of massive stellar 
populations, as a way to characterize the ratio of the ion- 
izing photon rate to the optical continuum. It is mainly 
sensitive to the global distribution of effective tempera- 
tures of the stars, although it also depends on the fraction 
of the ionizing flux absorbed by the nebular gas. Since this 
fraction is a priori unknown, we decided only to require 
the model to predict W(Hf3) greater than or equal to the 
observed value. 

The value of JU(H/3) used here is the ratio of 
the total nebular H/3 flux to the stellar H/3 contin- 
uum integrated over the cluster. The integrated 
values of these two quantities were obtained by 
dividing the values measured in the blue CAHA 
spectrum (cleared of the nebular continuum; see 
Section 13 .311 by the aperture throughputs of the 
optical slit for the nebular gas and for the stellar 
continuum (Section I3.4fl . This estimate of W(H(3) 
is affected by the error in the nebular-to-stellar ra- 
tio of the two aperture throughputs, evaluated to 
wlO%. We found W(R0) = 330 ± 30 A. 

4.3. UV lines 

The UV lines that for m in the atmospheres of mas - 
sive stars were studied bv lSekiguchi fc Anderson! l|l987a|) . 
They showed the correlation existing between the shapes 
and intensities of the SilV A1400 and CIV A1550 lines 
and the spectral type and luminosity class of early-type 
stars. The use of these two lines and, in some cases, the 



Hell A1640 line, can be a good complement to TU(H/3), 
for the determination of the evolutionary status and IMF 
of the cluster, due to the ir sensitivity to the luminos- 
ity class of the sta rs (e.g.. ISekiguchi fc Andersonlll987b[ 
iRohert et allll 993tlMa.s-Hesse fc Kuut^ll99llll99Qh . 

The IUE short-wavelength spectrum of NGC 588 has 
low resolution and signal-to-noise ratio, and was not used 
directly in the fit of the population parameters. Instead, 
it served to discuss and possibly reject some models de- 
termined with the other observables, by visual compari- 
son between model and observed spectra. Furthermore, we 
assumed this spectrum to be representative of the whole 
cluster, since the IUE slit covers nearly all the flux of the 
cluster (see Section EPjl . 

4.4. WR features 

The optical spectra clearly show the broad Hell A4686 
emission line produced by one or more WR stars. No corre- 
sponding CIV A5812 feature was found, showing that the 
detected WR stars are of WN type. The presence of such 
stars in the cluster fixes its age to 3.2-4.4 Myr if Z=0.004, 
and 2.8-4.5 Myr if Z=0.008, according to Starburst99. 
These age ranges are valid for m up = 120 Mq, and are 
narrower for lower mass cutoffs. 

Since the optical slit covers only a fraction of the clus- 
ter, our optical spectra do not necessarily include the sig- 
natures of all the kinds of stars present in the cluster. 
Consequently, we were unable to estimate the total num- 
ber of WN stars in the entire cluster. Therefore, to con- 
strain the models, we only used the fact that this number 
is of at least 1. Likewise, we ignored the non-detection of 
WC star in the optical spectrum, and did not process the 
number of WC stars predicted by the models. 

4.5. Model fitting and results 

We performed a chi-squarc-like fit of r, a and Z over 
W(H/3) and the number of WN stars 7V(WN), using the 
outputs of the Starburst99 code. A/"(WN), the total mass 
M to t and the number A/"(0) of O stars were multiplied by 
the coefficient needed to reproduce the luminosity of the 
cluster in filter F439W. The chi-square formula we used is 
the following: 

2 /max(330 -TU(H/3),0)\ 2 
X =\ 30 ) 

-2 In (l - e^Wj (i) 

The second term of the right part of this equations comes 
from the fact that the models described in this section 
predict numbers of WN stars, AA(WN), that are not in- 
tegers. Then, we treated A/"(WN) as the expectancy of a 
Poisson statistic process, and considered the likelihood to 
obtain at least one WN star with such a process, which is 
Pwn = 1 - exp(-A/"(WN)). 

The results of the fit are summarized in Table |3J 
The model associated to the minimum \ 2 is described 
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Fig. 3. Observed SED vs. SED constructed with the best- 
fit parameters. The observed SED is here supposed 
to be representative of the whole cluster, even 
though the slits (in particular the optical one) only 
sample a portion of the cluster. 





Z=0.004 


Observed 


Z=0.008 


W(H0) (A) 

7V(WN) 
reduced \ 2 


204 
0.17 
10 


330±30 
>1 


247 
0.03 
7 


t (Myr) 
a 

m up (Mq) 


3.5 (3.4-3.6) 
2.35 (2-3.1) 
110 (>90) 




2.8 (2.8-3.0) 
2.35 (2-3) 
120 (>110) 



Mtot (M ) (^3000) («3000) 
Af(O) («10) («12) 

Table 3. Results of the canonical analytical modeling of 
the cluster, for the two metallicities Z=0.004 and Z=0.008. 
Single values correspond to the best fits. Parenthesized 
ranges correspond to the approximate 3a ranges. 



by t =2.8 Myr, m up =120 M and Z=0.008. However, 
this model, with a reduced y 2 as high as 7, shows signif- 
icant discrepancies with the observations: both A/"(WN) 
and W(H/3) are well below the observed value. We can 
also note that due to the small value of jV(WN), the model 
Hell A4686 WR bump is negligible (see Fig. EJ - The bad 
fit of W(Hf3) when 7V(WN) is high enough can be at- 
tributed to the presence of blue supergiants (BSGs), de- 
fined here as having effective temperatures T e g < 3 x 10 4 
K, in the model cluster such that each of them, being alone 
in the nebula, would cause W(K0) < 100 A: according 
to the best-fit model, their number is jV(BSG) = 0.8 = 
0.07A/"(O). In a general way, such stars were predicted at 
any age with high enough values of A/"(WN), resulting in 
the absence of a good compromise between VF(H/3) and 
Af(WN); this is illustrated by Fig. H and 

So far, we were unable to determine the IMF slope 
more accurately than constraining it to the range 1.5- 
3. Now, we show how we can discard IMFs flatter than 
a r* 2. Fig. H3 shows the observed and model UV spectra 
in the case of the IMF slopes a — 1, 1.7, 2.35 and 3, and 
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Fig. 4. Evolution of W(R(3) (top), 7V(WN) (middle) and 
A/"(BSG) (bottom) with the cluster age, for Z=0.004, 
m up = 120 M Q and the IMF slopes a =2.35 (full lines), 1 
(short-dashed lines) and 3.3 (dotted lines). In the W(Rf3) 
panel, the full horizontal line shows the observed value, 
and the long-dashed one, the 3<r limit. In the 7V(WN) 
panel, the horizontal line and the arrow are here to re- 
mind that at least one WN star is present in the cluster. 



for the best-fit metallicity (Z=0.008), age (2.8 Myr) and 
upper mass limit (120 Mq). The synthetic UV spectra 
were created with the LMC/SMC library of Starburst99 
l|Leitherer et alJl200ll) . The depths of the SilV A1400 and 
CIV A1550 lines were found to decrease with increasing 
IMF slope, due to the conservation of their equivalent 
widths and increas e of their spec t ral wi dths. This result is 
in agreement with lRobert et all l)l993|) . Due to excessive 
depths of the lines, the cases a = 1 and 1.7 were rejected, 
while we considered the fits reasonably good for the other, 
steeper IMF slopes. 

5. Need to account for statistical fluctuations 

In the previous section, we obtained a model attempting 
to reproduce as well as possible W{R(3), the numbers of 
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Fig. 5. Same as Fig. [2 for Z=0.008 



WN stars and the UV stellar lines in the framework of 
a classical cluster analysis, that only takes into account 
the accuracy of the observational data. According to this 
model, the mass of the cluster is «3000 Mq. Meanwhile, 
assuming a meta llicity Z=0.004 and u s ing th e same evo- 
lutionary tracks, iMas-Hesse fc Kunthl l|l999|) obtained a 
cluster age r = 2.8, an IMF slope a = 1 and a total mass 
of 534 Mq , by the analysis of, mainly, the IUE-SWP and 
IUE-LWR spectra, but in absence of knowledge about the 
WR content of the cluster. 

In both works, the results have been obtained un- 
der the following assumptions: correct stellar evolution- 
ary tracks, instantaneous burst hypothesis and the im- 
plicit existence of a zero age main sequence (but see 
iTenorio-Tagle et al-lEoO^ and a number of stars in the 
cluster large enough to perfectly sample the IMF. However 
the estimated masses in both cases are about two orders of 
m agnitude below the 10 5 M^ mass for which, according 
to lCervino k Mas-Hessel l)l994l) . the whole IMF would be 
well filled. Hence, it is expected that sampling effects in 
the IMF will strongly affect any classical integrated anal- 
ysis of this cluster. 




1 .500 1 bbC 1 400 1 4bC 1 bOO 1 bb 
Wavelength (A) 




350 1400 1450 15 
Wavelength (A) 





dQ 1400 1450 1 500 1550 
Wavelength (A) 



(d) = 5 



1300 1350 1400 1450 15 
Wovelength (A) 



Fig. 6. Observed vs. synthetic UV spectra for different 
IMF slopes. The synthetic spectra were convolved by a 5 
A FWHM Gaussian curve to match IUE's resolution. 



Furthermore, the F439W luminosity of the cluster, 
with a value of -10.08 mag. is lower than the luminosity 
of the brightest s tar in the isochrone of the fi tted model, 
-11.42 mag., that ICerviho fc Luridianal l|2004h showed to 
be the lowest luminosity limit under which a cluster can- 
not be modeled with classical population synthesis. Even 
worse, this limit is brighter than the entire cluster by as 
much as 1.4 magnitude. This implies that the results of 
classical synthesis models regarding magnitudes and color s 
can be strongly biased llCervino fc Valls-Cabaiid! 120031) . 
Unfortunately, the actual theoretical status of statistical 
modeling of stellar clusters is not ready yet to solve such 
a situation. The only realistic assessments that can be ro- 
bustly obtained are: (i) due to the presence of at least one 
WN star, the age range is between 2.8 and 4.5 Myr (as- 
suming there is no issue in isochrone computations and 
that WR star formation is not due to the evolution of bi- 
nary systems), and (ii) the amount of gas transformed into 
stars at the onset of the burst does not exceed about 10 4 
Mq (for a Salpeter IMF and m\ ow = 1 Mq), for which the 
luminosity of the cluster would reach the lowest luminosity 
limit. 

This implies that the discussion about ages, IMF and 
discrepancies between the models and the observations in 
the analysis made with standard synthesis models hardly 
makes sense in our case. 

One of our main objectives is to infer a plausible shape 
of the ionizing continuum, which implies the use of a tech- 
nique different from the standard modeling as we show in 
the next section. 



6. Star-by-star characterization of the cluster 

Due to the small number of stars influencing the observ- 
ables and the ionizing flux of the NGC 588 cluster, and be- 
cause of the unsatisfactory results of the analysis detailed 
in Section^] we decided to study not only the cluster as a 
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whole, but also its individual stars. We achieved photom- 
etry of these stars, used isochrone curves to parameterize 
the cluster, attributed a model to each star, and synthe- 
sized spectra to be compared to the observations. 

6.1. Stellar photometry 

We performed photometric measurements of the stars on 
the five HST images, using the noao . digiphot . daophot 
package of IRAF. In each band, we ran daof ind and 
phot to detect candidate stars and measure their in- 
tensities in 6 pixel radius apertures. The brightest iso- 
lated stars served to compute analytical PSFs for all the 
bands. We then ran allstar to reject the artifacts of the 
daof ind routine and measure more accurately the inten- 
sities of the remaining objects. We selected again bright 
and isolated stars to carefully calculate the aperture cor- 
rections for the bands with the mkapf ile routine of the 
noao . digiphot .photcal package, and derive the stellar 
intensities integrated over the full PSF extents. The inten- 
sities were corrected for charge-transfer inefficiency and fil- 
ter contamination (see § 12. 4(1 . and converted into apparent 
magnitudes with the DN-to-flux keyword PHOTFLAM of 
the image headers,. The apparent magnitudes were con- 
verted into absolute (b ut reddened) ma gnitudes with the 
distance modulus 24.52 IjLee et all2002|) . The densest part 
of the cluster was not resolved by daof ind, but from vi- 
sual inspection, we found it to contain five clearly dis- 
tinct stars. Using the Levenberg-Marquardt's least-square 
method, and already knowing the PSF mathematical mod- 
els, we fitted the positions and fluxes of these five stars in 
each of the images; we detected no other star in the resid- 
uals of the fits, and retained our five-star detection. The 
final photometric lists were composed of 536 stars in band 
F547M, 698 in F439W, 353 in F336W, 574 in F170W and 
179 in F469N. 173 of them were common to bands F547M 
and F439W, 56, to F547M, F439W, F336W and F170W, 
and 20 to all bands. Table 01 lists the 56 stars common 
to F547M, F439W, F336W and F170W. In terms of abso- 
lute magnitudes, the 3<t detection limits were estimated to 
be approximately 0.8, —0.1, —0.8, —3.8 and —2.8 in filters 
F547M, F439W, F336W, F170W and F469N, respectively. 

The star observed with STIS is star #1 from Table H 
From the stars with available M^g, star #2 is the only 
one exhibiting a clear Hell excess, and was identified as 
the only or main WN star responsible for the Hell bump 
in the optical spectrum. 

6.2. Identification of the star observed with STIS 

To identify the nature of the star observed with STIS, 
we first compared it s spectrum wit h the ones of the 
LMC/SMC library of lLeitherer et alJ ll200lh . paying at- 
tention to the SilV A1400, CIV A1550 and, when available, 
Hell A1640 lines. The four acceptable matches (Fig. [TJ) in- 
dicated the star to be of class BOII, BOI, late WN (WNL) 
or late WC (WCL). We found a significant Hell A4686 








1300 '400 1500 '600 1700 '300 1400 '500 

Wavelength (A; Wavelength (A) 



Fig. 7. STIS vs. LMC/SMC spectra: BOII, BOI, WNL and 
WCL, respectively The 1600-1700 A region was unavail- 
able for the BOI star. 
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Fig. 8. Optical spectra near the WR Hell A4686 bump, 
extracted from the blue CAHA spectrum around stars #1 
and #2 of Table H 



bump in the CAHA slit zone situated within the PSF ex- 
tent of this quite isolated star (see Fig. [S] for this star and 
for star #2 of Table 0J. The intensity of this bump be- 
ing compatible with the (quite inaccurate) value derived 
from the photometric measurements, we deduced that the 
star is a WR. In absence of any optical carbon feature, 
in contrast with the well-defined Hell A4686 bump, we 
concluded that it is a WNL star. 

Considering the presence of such a star in the cluster, 
we inferred the age range of the latter for both metal- 
licities Z=0.004 and Z=0.008: 3.2-3.7 and 2.8-4.5 Myr, 
respectively. 
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Star # 


R.A. 

01h32m 


Dec. 
+30° 


M547 




M439 


M 33 6 




M170 


M 4 69 


Eb-v 


mini 

(M ) 


1 


45.23s 


38'58.4" 


-6.63±0 


04 


-7.42±0.03 


-8.14±0 


04 


-9.34±0.04 


-7.26±0.06 


0.11±0.03 


55.3 


2 


45.56s 


38' 54.4" 


-6.03±0 


04 


-6.91±0.04 


-7.74±0 


03 


-9.32±0.04 


-6.68±0.08 


0.00±0.03 


55.5 


3 


45.50s 


39'07.0" 


-6.12±0 


04 


-6.91±0.03 


-7.55±0 


04 


-8.57±0.05 


-7.09±0.06 


0.11±0.02 


43.8 


4 


45.31s 


39'05.7" 


-5.60±0 


05 


-6.45±0.04 


-7.08±0 


04 


-8.24±0.07 


-6.19±0.09 


0.11±0.03 


40.1 


5 


46.81s 


39'06.9" 


-5.34±0 


07 


-6.22±0.05 


-6.89±0 


06 


-8.26±0.08 


-5.91±0.14 


0.08±0.04 


38.1 


6 


45.37s 


38'53.4" 


-5.34±0 


07 


-6.19±0.07 


-6.89±0 


06 


-8.33±0.06 


-5.57±0.13 


0.08±0.04 


37.1 


7 


45.58s 


38'55.5" 


-4.86±0 


05 


-5.90±0.05 


-6.55±0 


05 


-8.29±0.06 


-5.02±0.08 


0.02±0.04 


32.1 


8 


45.37s 


39'06.4" 


-4.80±0 


07 


-5.59±0.06 


-6.45±0 


06 


-7.73±0.10 


-5.36±0.15 


0.11±0.02 a 


33.1 


9 


45.53s 


38'54.9" 


-4.64±0 


07 


-5.49±0.06 


-6.35±0 


05 


-7.96±0.07 


-5.04±0.15 


0.07±0.01 a 


32.1 


10 


45.60s 


38'55.6" 


-4.56±0 


06 


-5.31±0.06 


-6.17±0 


06 


-7.46±0.09 


-4.40±0.08 


0.12±0.01 a 


30.1 


11 


45.59s 


38'55.4" 


-4.27±0 


07 


-5.26±0.06 


-5.78±0 


09 


-7.46±0.09 


-4.85±0.09 


0.09±0.01 a 


29.1 


12 


45.47s 


38'55.7" 


-4.42±0 


06 


-5.19±0.05 


-6.02±0 


06 


-7.51±0.08 


-4.94±0.15 


0.10±0.01 a 


29.1 


13 


45.55s 


38'55.3" 


-4.28±0 


07 


-5.03±0.07 


-5.93±0 


07 


-7.39±0.09 


-4.77±0.17 


0.10±0.01 a 


27.1 


14 


45.59s 


38'55.7" 


-4.02±0 


09 


-4.90±0.07 


-5.39±0 


09 


-6.93±0.13 


-4.79±0.09 


0.12±0.02 a 


26.1 


15 


46.07s 


39'02.5" 


-4.10±0 


06 


-4.79±0.06 


-5.51±0 


07 


-6.54±0.14 




0.18±0.02 


27.1 


16 


45.60s 


38'55.4" 


-3.87±0 


08 


-4.75±0.07 


-5.37±0 


08 


-6.93±0.11 


-4.45±0.09 


0.10±0.02° 


24.4 


17 


45.80s 


39'02.0" 


-3.71±0 


07 


-4.50±0.07 


-5.32±0 


07 


-6.63±0.13 


-3.87±0.27 


0.11±0.02 a 


21.9 


18 


44.61s 


38'52.6" 


-3.35±0 


11 


-4.28±0.08 


-4.88±0 


09 


-6.50±0.12 


-3.04±0.50 


0.08±0.02 a 


20.1 


19 


45.42s 


39'08.5" 


-3.36±0 


08 


-4.22±0.08 


-5.00±0 


08 


-6.16±0.17 


-3.33±0.42 


0.12±0.02 a 


19.7 


20 


45.05s 


38'55.2" 


-3.28±0 


08 


-4.19±0.07 


-5.01±0 


09 


-6.66±0.11 


-3.60±0.36 


0.05±0.02 a 


19.5 


21 


45.54s 


38'55.9" 


-3.31±0 


11 


-4.14±0.08 


-4.86±0 


10 


-6.44±0.14 




0.08±0.02 a 


19.1 


22 


45.20s 


38'59.3" 


-3.19±0 


13 


-4.07±0.10 


-4.72±0 


14 


-6.08±0.19 




0.11±0.03° 


18.6 


23 


45.10s 


38'58.9" 


-3.05±0 


08 


-3.99±0.08 


-4.71±0 


10 


-6.13±0.13 


-3.50±0.35 


0.08±0.02 a 


18.1 


24 


45.67s 


38'57.4" 


-3.09±0 


09 


-3.97±0.09 


-4.55±0 


12 


-6.08±0.14 




0.09±0.02 a 


17.9 


25 


45.16s 


39'04.2" 


-3.05±0 


11 


-3.87±0.11 


-4.15±0 


15 


-5.41±0.20 




0.18±0.03 a 


17.2 


26 


45.24s 


38'56.2" 


-2.94±0 


09 


-3.81±0.09 


-4.37±0 


12 


-5.94±0.15 




0.09±0.02 a 


16.9 


27 


45.36s 


38'54.6" 


-2.86±0 


09 


-3.76±0.08 


-4.42±0 


11 


-5.87±0.16 




0.09±0.03 a 


16.6 


28 


45.21s 


38'55.8" 


-2.85±0 


11 


-3.72±0.09 


-4.28±0 


14 


-5.65±0.21 




0.12±0.03 a 


16.3 


29 


46.08s 


39'01.7" 


-2.93±0 


10 


-3.71±0.10 


-4.03±0 


13 


-4.71±0.32 




0.25±0.04 


20.0 


30 


45.20s 


38' 54.2" 


-3.07±0 


18 


-3.62±0.10 


-4.12±0 


14 


-5.74±0.19 




0.13±0.04 a 


15.6 


31 


45.62s 


38'58.6" 


-2.75±0 


10 


-3.57±0.09 


-4.47±0 


11 


-5.91±0.20 




0.06±0.03 a 


15.5 


32 


45.46s 


38'54.5" 


-2.64±0 


10 


-3.51±0.09 


-4.26±0 


12 


-5.69±0.19 




0.08±0.03 a 


15.1 


33 


45.85s 


39'07.1" 


-2.52±0 


11 


-3.34±0.10 


-4.07±0 


12 


-5.36±0.24 




0.10±0.04 a 


14.2 


34 


45.33s 


38'56.7" 


-2.39±0 


10 


-3.33±0.10 


-4.10±0 


12 


-5.73±0.17 




0.03±0.03 a 


14.3 


35 


45.80s 


39'11.5" 


-2.34±0 


11 


-3.30±0.09 


-4.01±0 


13 


-5.67±0.19 




0.03±0.03 a 


14.1 


36 


45.72s 


39'10.6" 


-2.33±0 


13 


-3.27±0.11 


-3.89±0 


13 


-5.34±0.23 




0.08±0.04 a 


13.9 


37 


45.51s 


38'55.9" 


-2.42±0 


10 


-3.22±0.09 


-3.70±0 


14 


~5.26±0.24 




0.10±0.04 a 


13.6 


38 


46.15s 


38'56.7" 


-2.32±0 


10 


-3.14±0.11 


-3.74±0 


13 


-4.85±0.28 




0.14±0.04 a 


13.2 


39 


45.22s 


38'54.7" 


-2.39±0 


16 


-3.09±0.14 


-3.77±0 


16 


-4.89±0.29 




0.14±0.05 a 


12.9 


40 


45.12s 


38'55.8" 


-2.40±0 


13 


-3.04±0.12 


-3.84±0 


15 


-5.12±0.24 




0.11±0.04 a 


12.7 


41 


45.62s 


38' 54.2" 


-2.27±0 


11 


-3.00±0.11 


-3.69±0 


13 


-5.07±0.25 




0.10±0.04 a 


12.6 


42 


45.62s 


38'57.8" 


-2.06±0 


13 


-2.99±0.12 


-3.67±0 


17 


-5.17±0.25 




0.05±0.04 a 


12.6 


43 


45.61s 


38'55.1" 


-2.39±0 


19 


-2.97±0.17 


-3.60±0 


19 


-5.06±0.30 




0.11±0.05 a 


12.4 


44 


45.33s 


38'55.8" 


-2.16±0 


13 


-2.95±0.12 


-3.70±0 


13 


-5.37±0.26 




0.03±0.04 a 


12.4 


45 


44.74s 


38'55.8" 


-1.78±0 


15 


-2.94±0.13 


-3.07±0 


20 


-5.11±0.29 




0.01±0.05 a 


12.4 


46 


45.33s 


38'53.0" 


-2.23±0 


12 


-2.91±0.12 


-3.66±0 


16 


-4.92±0.28 




0.11±0.05 a 


12.1 


47 


45.27s 


38'58.6" 


-1.98±0 


14 


-2.89±0.13 


-3.56±0 


17 


-5.11±0.25 




0.04±0.04 a 


12.1 


48 


45.58s 


38'55.8" 


-2.32±0 


20 


-2.88±0.19 


-3.94±0 


13 


-5.25±0.26 




0.07±0.05 Q 


12.1 


49 


45.01s 


38'50.9" 


-1.90±0 


15 


-2.84±0.12 


-3.54±0 


18 


-5.34±0.28 




0.00±0.05 a 


12.0 


50 


45.39s 


38'55.6" 


-1.91±0 


15 


-2.83±0.14 


-3.50±0 


16 


-4.92±0.30 




0.05±0.05 a 


11.8 


51 


45.23s 


38'54.9" 


-2.04±0 


14 


-2.78±0.15 


-3.37±0 


15 


-4.34±0.44 




0.16±0.07 a 


11.5 


52 


45.18s 


38'53.0" 


-1.81±0 


15 


-2.70±0.13 


-3.26±0 


16 


-5.06±0.26 




0.01±0.05 a 


11.4 


53 


44.86s 


38'57.5" 


-1.60±0 


17 


-2.59±0.14 


-3.01±0 


23 


-4.79±0.31 




0.01±0.06 a 


10.9 


54 


45.67s 


38' 54. 7" 


-1.67±0 


15 


-2.49±0.17 


-3.18±0 


19 


-4.77±0.29 




0.02±0.05 a 


10.5 


55 


45.66s 


38'53.2" 


-1.36±0 


18 


-2.27±0.17 


-2.94±0 


19 


-4.48±0.42 




0.00±0.07 a 


9.6 


56 


44.96s 


38'55.7" 


-1.30±0 


19 


-2.09±0.19 


-2.84±0 


20 


-4.43±0.48 




0.00±0.08 a 


9.0 



Table 4. Positions, reddened absolute magnitudes, Eb-v coefficients, and initial masses of the brightest stars of NGC 
588. The extinction is discussed in § 16.3.21 and the initial masses, in § 16.41 "Values eventually set to 0.11 
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6.3. Comparison with isoch rones 

6.3.1. Construction of the isochrones 

We constructed isochrone curves calling the same evolu- 
tionary tracks and model atmospheres as in Section 0] 
The effective temperatures and bolometric luminosities 
of these isochrones were computed with the original 
Starburst99 code. The corresponding magnitudes, how- 
ever, were derived from them by means of an interpola- 
tion method different from the one of Starburst99. The 
latter mainly consists in nearest-neighbor selection of the 
spectra in the (logT e ff, logg) plane, T c s being the effec- 
tive temperature, and g, the surface gravity. This induces 
step-like discontinuities in plots such as color-magnitude 
diagrams. This is why we exploited an alternative interpo- 
lation method, explicated in the Appendix, that consists, 
at each wavelength, of applying a bilinear interpolation of 
the logarithm of the flux in the (log g,y) plane, y being a 
variable that depends both on the effective temperature 
and on the wavelength. 

Depending on the stellar parameters, various kinds of 
model atmospheres, listed in Section 0] were exploited. 
This implies an artificial jump in the isochrone dia grams 
at the locus of transition fromlLeieune et alJ l|l997|) mod- 
els to the one of lSmith et alJ l|2002() . Fortunately, this tran- 
sition occurs in a zone where it is smaller than the error 
bars, and consequently it has a negligible impact. 

6.3.2. Reddening 

We first investigated the extinction law, since a priori it 
may be different from the one inferred from the canon- 
ical study of the cluster (cf. § 14. If) . We considered sev- 
eral isochrones in the age range of the cluster, and for 
each of them and each of the three tested extinction laws 
(Galactic, LMC and SMC), we proceeded as follows. We 
first estimated the color excess of each star in Table 0] 
by shifting the measured points towards the isochrone 
curve in the (M4 39 ,M 336 — M547) diagram according to 
the tested law, assuming stars #1 and #2 to be WN stars, 
and all the other ones to be on the main sequence. Then, 
we dereddened each star in the (A/439, A/170 — -M547) di- 
agram, and compared the resulting points with the theo- 
retical isochrone. Independently of the reference isochrone 
(i.e., in practice, the age), the SMC law was favored, as 
the systematic discrepancy between the isochrone and the 
locus of the stars, in particular the brightest ones, was 
small for the SMC law and significantly larger for the 
other two laws. Numerically, the discrepancy between the 
observational points in the (A/439, A/170 — -M547) diagram 
dereddened with the use of the (M439, M336 — M547) one 
and the theoretical isochrone is characterized by a reduced 
X 2 ranging in the approximate intervals 2.3-3.3 for the 
Galactic law, 1.8-2.0 for the LMC one and 0.6-0.8 for the 
SMC one, obtained with 56 stars. The fact that this % 2 
is close to 1 with the SMC law made us confident in the 
latter, and is also an indication that none of the 56 con- 




-2.0 -1.B -1.6 -1.4 
M336-M547 





Fig. 9. Illustration of the dereddening procedure. The 
isochrone shown here was established for an age of 3.5 Myr 
and for Z=0.008. Panel (a): observed ( AI439 , M 336 - A/547) 
diagram. Other panels: (AT439, A/170 — A/547) diagram 
dereddened with the (A/439, A/ 336 — A/ 547 ) one and the 
Galactic (b), LMC (c) and SMC (d) laws. The dotted 
lines are the WR branches of the isochrone curves. The 
two WR stars are signalized by open circles. The arrows 
indicate the shift that the observational points undergo by 
if dereddened for Eb-v = 0.1. In the diagrams of panels 
(b), (c) and (d), the error bars account for uncertainties 
of E B -v- 




Fig. 10. Same as Fig. El for the isochrone of 4.5 Myr. 
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20 - 




X=(Eb-v-0.11)/u{E B -v) 

Fig. 11. Histogram of the dispersion of the stellar color 
excesses Eb-v around the nebular one 0.11. The ab- 
scissa axis variable is the discrepancy term \ = (Eb-v ~ 
0. 11) / a (Eb-v), where a (Eb-v) is the uncertainty of 
the considered measurement of Eb-v- The histogram in- 
cludes the 173 stars observed simultaneously through fil- 
ters F547W and F439W. The overplotted curve is the the- 
oretical expected distribution in the case of a unique real 
value Eb-v =0.11 and of well-determined measurement 
uncertainties. 



sidered stars is significantly affected by blending with an 
unresolved companion. The procedure of selection of an 
extinction law is illustrated in Fig. [§] for an age of 3.5 
Myr, and in Fig.UHlfor 4.5 Myr. 

Once the extinction law chosen, we re-considered the 
individual color excess Eb-v of the stars. For the stars 
used in the selection of the extinction law, we started from 
the values already estimated in the (M439,M336 — M547) 
diagram. For the stars that, once dereddened, were found 
to be brighter than M439 = —6, we conserved the values of 
Eb-v found in this diagram and the corresponding uncer- 
tainties, the latter including the (small) scatters of values 
found with the different possible isochrones. However, for 
the stars dimmer than M439 = —6, we used the constancy 
of the isochrones in the (M439, M170 — M547) diagram to 
estimate Eb-v with better accuracy. Indeed, the uncer- 
tainty in Eb-v determined in one of the (M439, Mx — 
M547) diagrams (X = 170 or X = 336) is roughly the ra- 
tio a(M x -M 5i7 )/( fx~hii), where a(M x -M 5i7 ) is the 
uncertainty in the Mx — M547 color, and in our case, this 
ratio in Eb-v is smaller for X = 170 than for X = 336. 
With a few exceptions, we found all the color excess to 
be compatible with the nebular value, Eb-v = 0.11, and 
when this was the case and the dereddened F439W mag- 
nitude was dimmer than —6, we set them to this value. 
An important departure from Eb-v = 0.11 is the case of 
star #2: E B -v = 0.00 ± 0.03. All stars with unavailable 
F170W and F336W magnitudes were dereddened with 
Eb-v = 0.11, compatible with the (M439,M439 — M547) 
theoretical curve. Fig.EJtestifies to the correctness of set- 




-5.0 -4.5 -4.0 -3.5 -3.0 -2.5 -5.C -4.5 -4.0 -3.5 -3.0 -2.5 

Ml70-Mt,47 Mi 7C -Mt,47 



Fig. 12. Observed vs. isochrone (M4 3 g,Mi7o — M547) di- 
agram for Z=0.004. The two WR stars are marked with 
open circles. The WR branches of the isochrones are shown 
as dashed lines. 




-4.5 -4.0 -3.5 -3.0 -2.5 -5.0 -4.5 -4.0 -3.5 -3.0 -2.5 
M17Q-M547 M17C-M547 




M170-M547 M17C-M547 
Fig. 13. Same as Fig. C3 for Z=0.008. 

ting almost all the 173 stellar color excesses to the nebular 
one. 

In all what follows, the color-magnitude diagrams 
(CMDs) refer to the dereddened stellar magnitudes. 

6.3.3. Fit of the age and of the metallicity 

The age and metallicity of the cluster were constrained 
by simultaneously fitting the isochrone to the dereddened 
observational points in the (M43g, M\ 70 — M547) diagram, 
and by comparing synthetic spectra and derived properties 
to the observations. For a given isochrone, each star was 
identified to the nearest point of the isochrone in terms 
of chi-square (i.e., accounting for the error bars) in the 
(M439, M170 — M547) diagram or, when the F170W magni- 
tude was unavailable, in the (M43g, M439 — -M547) diagram, 
thus attributing the most appropriate model spectrum to 
the considered star. If the cluster spectrum was synthe- 



L. Jamet et al.: The ionizing cluster of NGC 588 



13 




-4 - 



Fig. 14. Final adopted (M439, Mi 70 — M547) diagram. 



HR diagram 



6 - 



5 - 



A - 



3 - 



5.0 




4.4 4.2 

log T.„ (K) 



Fig. 15. Final adopted HR diagram for all the detected 
stars. The circles show the selected individual stellar 
model. The dashed part of the diagram is the WR branch. 
The two dotted curves represent the HR diagram at 3.0 
Myr (upper curve) and 5.5 Myr (lower curve) 



sized to simulate what is seen through a slit (CAHA or 
IUE), the spectrum selected for a given star was multi- 
plied by its corresponding aperture throughput (known 
from its position with respect to the slit and from the 
angular PSF or seeing of the considered data) and a func- 
tion representing its possible differential extinction (with 



respect to the reference value Eb-v = 0.11) along the 
processed wavelength range, and added to the total clus- 
ter spectrum to compute. If the total cluster spectrum was 
to be synthesized, then the individual stellar spectra were 
summed without being previously modified. 

We computed the isochrone curves for the different 
ages multiple of 0.1 Myr belonging to the ranges con- 
strained by the presence of a WNL star. Here, we show 
the curves obtained for the ages 3.2 and 3.7 Myr in the 
case Z=0.004 (Fig. 02), and 3.0, 3.5, 4.0 and 4.5 Myr for 
Z=0.008 (Fig. EJ>- At both metallicities, the isochrone 
high-luminosity main-sequence branch was found to move, 
with increasing age, from regions of low values of the 
M170 — M547 color towards the somewhat "red" obser- 
vational points. 

The synthesized spectra were exploited as follows. For 
a given tested isochrone, the spectrum simulated for 
the optical slit was used to extract the model ratio 
R B = F A (3630)/F A (3780) related to the Balmer jump, 
which is a diagnostic of the effective temperature of a clus- 
ter dominated by hot stars. Indeed, in the spectrum of a 
hot star, the intensity of the Balmer jump decreases with 
increasing effective temperature, and usually serves to de- 
termine the subtypes of B stars, whereas it is negligible 
in the spectra of O stars. The observed value, whose un- 
certainty is dominated by the local fluctuating residues in 
photometric calibration, is Rb = 1.12 ± 0.08. The con- 
structed UV spectrum was visually compared to the IUE 
observed one. Finally, we derived the predicted W(Hf3) 
from the total cluster spectrum, assuming again that the 
nebula absorbs all the ionizing photons, and thus request- 
ing the predicted W(H(3) to be greater than or equal to 
the observed one. 

Mathematically, we used the following \ 2 estimator to 
constrain the age and the metallicity of the cluster: 



max(330 -IU(H/3),0) 
30 



R B - 1-12 
0l)8 



-v? 

AlSOC 



(2) 



where Xisoc characterizes, for a given isochrone, the dis- 
crepancy between the latter and the loci of the brightest 
four MS stars in the (M43g,Mi7o — M547) diagram. The 
numerical results are summarized in Table |3J In prac- 
tice, the synthetic Balmer jump was constant, and the 
fit was constrained by W(Kf3) and xf soc . For Z=0.004, 
even at 3.7 Myr, the theoretical isochrone remained to- 
ward values of Mi 70 — M547 lower ( "bluer" ) than the ob- 
servational points of the brightest stars, though the fit 
was acceptable. Meanwhile, for Z=0.008, the fit of the 
isochrone was very satisfactory for ages around 4.0 Myr, 
with W(R/3) also being well predicted. We retained the 
model with Z=0.008 and r = 4.2 Myr as the best-fit 
one. The (M43g,Mi7o — M547) diagram for this model 
is shown in Fig. 1141 and the corresponding Hertzsprung- 
Russel (H R) diagram, in F ig. 115] According to the classi- 
fication of ISchmidt- Kaler (1982), 20 stars were classified 
as O-type ones. 
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Z=0.004 Observed Z=0.008 

xL~c 5^ 2~3 

W(R/3) (A) 602 330±30 342 

R B 1T0 1.12T0.08 1.10 

age (Myr) 3.7 (3.7-3.7) 4.2 (3.6-4.4) 

Table 5. Numerical results of the star-by-star analysis of 
the cluster. The parenthesized age ranges indicate the 90% 
confidence limits. 

Quantity L(HeII A4686) NIII A4640/HeII A4686 

(10 35 erg s" 1 ) 

Observed 3.0±0.2 0.17±0.01 

Model 13±9 0.2±0.1 

Table 6. WR line strengths, observed/predicted in the 
dereddened CAHA spectrum. 

The SED of the best-fit model is shown in Fig. ^] 
We considered it to satisfactorily reproduce the observa- 
tions, except for the strength of the optical WR bumps. 
However, the model atmospheres used for WR stars here 
are not intended to compute reliably these bumps, and 
we did not pay attention to the ir strength. Instead , we es - 
timated them from Table 1 of ISchaerer fc Vaccal l|l99§(l . 
knowing that star #2 was finally classified as a WNL 
from its hydrogen surface abundance. The results, satis- 
fying but somewhat inaccurate (due in particular to the 
scatter of the model bump intensities), are summarized 
in Table |H| The synthetic UV spectrum, superimposed to 
the observed one in Fig. El wa s found to satisfactorily 
reproduce the CIV A1550 line. It also contains an unde- 
tected SilV A1400 feature. However, the UV line spectrum 
is dominated by a small number of massive stars, that are 
not necessarily well represented by the spectra LMC/SMC 
library, in particular the WNL ones. In consequence, we 
did not consider the SilV A1400 discrepancy as a critical 
one. 

6.4. Mean IMF 

From the initial mass associated with each star, we com- 
puted the mean power-law IMF of the cluster. More specif- 
ically, we constructed a histogram of the logarithm of the 
initial masses, from log(m) = 0.8 to log(m) = 1.76 (m = 
6.3 to m = 58 Mq), with a bin of 0.03 (equivalent to a fac- 
tor of 1.07 between two successive bins), and fit it with an 
exponential law, assuming Poissonian noise in each bin, 
knowing that a power-law IMF dN/dm oc mT a can be 
translated into the law dN/dlog{m) oc lot 1 -") 1 ^™). The 
lower limit of the initial mass range was chosen in order 
to ensure that the stars of the cluster belonging to a given 
bin were all detected, and the upper limit corresponds to 
the highest possible initial mass at the age and metal- 
licity of the cluster. The result is shown in Fig. ^] We 
found an IMF slope a — 2.37 ± 0.16, and retained the 
compatible Salpctcr slope, resulting in an inferred IMF 
dN/dm = (2680 ± 210) to" 2 ' 35 . The corresponding total 
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2000 4000 6000 8000 

Wavelength (A) 



Fig. 16. Optical and UV SEDs synthesized for Z=0.008 
and r = 4.2 Myr with the individual stellar models. For 
better visibility and comparison with Fig.0 the luminosi- 
ties of the observed optical and UV parts were rescaled 
with multiplicative factors, in order to fit the integrated 
luminosities in band F439W and F170W, respectively. The 
synthetic spectra were normalized using the same factors. 
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Fig. 17. UV spectrum synthesized for Z=0.008 and r = 
4.2 Myr with the individual stellar models. 

initial mass integrated in the range 1 < m < 58 Mq is 
Mtot = 5800 ± 500 M . 

7. Discussion 

The most massive stars of an OB association are also the 
most influential ones, though the least numerous, on the 
spectrum of this cluster. In moderately massive clusters, 
their small number is subject to significant fluctuations 
around the mean IMF, causing large variations in the 
observable (near UV and optical) and Lyman continuum 
spectral ranges. For instance, in Section l4~5l we saw the 
significance of the presence of only 0.8 BSG in the analyt- 
ical model of NGC 588 cluster upon the overall properties 
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100.0 



L 




0.1 I , , , I , , , I , , , I , , , III, II, 

0.8 1.0 1.2 1,4 1.6 

log(m) (M Q ) 

Fig. 18. IMF histogram and fit. Full line: best fit. Dashed 
line: fit for the Salpeter slope. The error bars are equiva- 
lent to la Gaussian ones, in terms of likelihood, and must 
be interpreted this way: for a model that goes through 
the extremity of an error bar, the observed value whose 
error bar is being considered has a likelihood to occur of 
exp(— 1/2) times the maximum likelihood of the model at 
this abscissa. 

derived from this model. We now discuss the effects of the 
real sampling of the IMF in various mass ranges associ- 
ated to different kinds of massive stars, and more gener- 
ally in the mass range covered by the most luminous stars. 
Then, we present the results of a simulation of the effect of 
discrete random IMF sampling, in the whole initial mass 
range, over W(H/3). 

7.1. Mean and observed numbers of different kinds of 
massive stars 

From the mean IMF computed in Section l6.4l we derived 
the mean expected numbers of several kinds of massive 
stars susceptible to change significantly the spectral prop- 
erties of the cluster. Here we present these results and 
discuss them. 

7.1.1. Blue supergiants 

According to our stellar models, at 4.2 Myr and for 
Z=0.008, BSGs (as defined in § 14.5(1 are stars with ini- 
tial masses ranging from 44 to 55 Mq. The number ex- 
pectancy, derived from the integration, between these two 
limits, of the mean power-law IMF computed in § 16.41 
is TV(BSG) = 3.0. The associated Poisson probability to 
observe no BSG, as is our case, is exp(A/"(BSG)) = 5%. 
We added an artificial continuous population composed 
by this kind of stars, weighted by the Salpeter IMF slope, 
to the actual stellar content of the cluster. Independently 
of the optical absolute luminosity - which would be au- 
tomatically fit in a classical modeling of the cluster - the 
most striking change in the computed spectral properties 



is, as expected, the large drop of W(H(3): in presence of 
the BSG population, the latter quantity would be only 
133 A, instead of 342 A. 

7.1.2. Wolf-Rayet stars 

We computed the expected number of the different WR 
stars (WNL, WNE, WCL, WCE, WO), derived from 
the mean IMF of | we found A^(WNL) = 0.15, 
AA(WNE) = 0.02, JV(WCL) = 0.02, TV(WCE) = 0.01 
and A/"(WO) = 0.29, according to the star classifi- 
cation of Starburst99. The resulting Poisson likelihood 
of the observed WR content of the cluster is 0.7%, a 
small value that is however much more satisfactory than 
the 0.04% likelihood derived from the analytical model. 
Furthermore, we simulated a spectrum obtained by re- 
placing the two WR stars of the cluster with an analyt- 
ical, complete population of 0.51 WN star following the 
Salpeter IMF. This population was expected to produce 
less luminous but possibly harder or softer radiation than 
the two actually detected WR, as the WR branch of the 
HR diagram (Fig. 115(1 spans a significant range of effective 
temperatures. We found that whereas W(H(3) is smaller 
with the modified stellar population, the latter produces 
harder radiation than the "true" one in the Lyman con- 
tinuum range, especially below the threshold of ionization 
of He + , as summarized in Table 

7.2. On the importance of the few most massive stars 

In § 17.1. II and !7. 1.21 we saw the significance of fluctuations 
in the population of some characteristic kinds of massive 
stars upon the spectral properties of a cluster as moder- 
ately massive as the one of NGC 588. This issue can be 
generalized to all sorts of very massive stars, that can radi- 
ate intensely in the optical and/or ionizing spectral range 
and whose numbers are subject to the largest relative sta- 
tistical fluctuations. 

In Fig. E| are plotted the SEDs of the whole clus- 
ter, of the 6 brightest stars (in filter F439W) and of the 
remaining 167 stars, as well as the fractional contribu- 
tion of the 6 brightest stars to the total SED. It is ev- 
ident that these 6 stars are very influential on the to- 
tal spectrum, whether at observable wavelengths (where 
they are responsible for about half the total flux) or in 
the high-energy range, where they generate approximately 
two thirds of the hydrogen-ionizing photons, even though 
they constitute only about 1/30 of the whole set of de- 
tected stars. These 6 stars, which also are the most massive 
ones, are specifically the stars situated in the temperature- 
luminosity zone that varies significantly within the first 
few millions of years of the cluster, as one can see in 
Fig. US Fig. H3 shows the (Af 439 , log Q(H )) diagram of 
the detected stars. The most striking feature is the quasi- 
vacuousness of the curve in the region of the 6 dominating 
stars, in particular the BSG branch that drops towards 
negligible Lyman continuum at high band F439W lumi- 
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BSG content 


real 


IMF 


real 


IMF" 


2.8 Myr a ' b 


WR content 


real 


real 


IMF 


IMF" 


2.8 Myr"' 6 


MS content 


real 


real 


real 


IMF" 


2.8 Myr a - b 


W(Uf3) (A) 


342 


133 


284 


80 


247 


Q(H°) (10 50 s" 1 ) 


2.14 


2.16 


1.39 


0.54 


1.60 


Q(He°)/Q(H°) 


0.13 


0.12 


0.14 


0.11 


0.14 


10 s Q(He+)/Q(H°) 


9.3 


8.0 


2.8xl0 4 


8.8 


130 


Q(0+)/Q(H°) 


0.020 


0.020 


0.024 


0.011 


0.017 



Table 7. Some spectral properties derived from the observed stellar content (2nd column) and from the one with BSG, 
WR and whole population replaced with the one predicted by the fitted power-law IMF (3rd, 4th and 5th columns, 
respectively). Q(X) designates the photon rate able to ionize the ion "X" population of a nebula. a The 5th and 6th 
columns correspond to analytical models, whose total mass here reproduce the observed F439W luminosity. fc Best-fit 
model of the canonical analysis. 
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Fig. 19. Top: SEDs of the whole cluster (full line), of the 
6 brightest stars (dashed line) and of the remaining 167 
stars (dotted line). Bottom: ratio of the brightest 6 stars 
SED (full line) and of the remaining 167 star SED (dotted 
line) to the total SED. 

nositics. If the IMF of the cluster were complete, as as- 
sumed in a classical model of stellar population, then the 
main spectral diagnostic, W(H{3), would be significantly 
lower than what we observe at the age and metallicity we 
derived from the star-by-star approach. This is illustrated 




log Q(H°) (s- 1 ) 

Fig. 20. (log Q(H°), M439) diagram of the detected stars 
along the the theoretical isochrone. The WR branch is the 
dashed part of the curve. 

in Fig. [21 where the SED found with the classical ap- 
proach is also shown. In this figure, we can appraise the 
significance of the drop of the Lyman-to-optical contin- 
uum ratio, as well as the increase of the magnitude of the 
Balmer jump and the decrease of the global slope of the 
SED at observable wavelengths, when passing from the 
star-by-star SED to the one of the analytical model. 

7.3. Approximate expected uncertainties and 
Monte-Carlo simulation 

In Section 17711 and [7"21 we saw the significant consequences 
of the deviation of the actual (discrete and statistically 
fluctuating) IMF of NGC 588 from the mean and contin- 
uously sampled analytical one. In this cluster of ~6000 
Mq (integrated in the interval of individual initial stellar 
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Fig. 21. Star-by-star constructed SED (full lines), and 
SEDs of the analytical populations, for Z=0.008 and 
m up = 120 M , at 2.8 Myr (dashed line) and 4.2 Myr 
(dotted line). 
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Fig. 22. Histogram of logW(H[3) resulting from the 10 6 
Monte-Carlo discrete IMF realizations. 

masses 1-58 Mq), the stars massive enough to rule its 
spectral properties, whether observable (W(H{3), Balmer 
jump, etc.) or not (e.g., hardness of the Lyman contin- 
uum), are only 6, which can easily explain the encountered 
issues. 

In order to estimate quantitatively the uncertainties 
related to the fluctuations of massive stellar populations 
of same age, metallicity and IMF as the ones of NGC 588, 
we performed a simple Monte-Carlo simulation where the 
mean IMF of § 16.41 was divided into 0.01 Mq wide mass 
bins, in which 1977 stars (1977 being the number of stars 
integrated in the mean IMF) were randomly distributed, 
according to the probability of each bin to be filled with a 
given star, this probability being given by the mean ana- 
lytical IMF. We performed 10 6 such Monte-Carlo realiza- 
tions. This model assumes the stars of the cluster to have 
random masses, uncorrelated to each other (a simplifica- 
tion that may be corrected when the actual properties of 



star formation in clusters are better known), and follow- 
ing a likelihood function given by the mean IMF. For each 
random distribution, we calculated the resulting W{Hj3)-, 
we found log(lF(H/3)) to follow a roughly Gaussian dis- 
tribution with stand ard deviation 0.24 d ex (as also ob- 
tained analytically bv lCervifio et alJ2002(l . The histogram 
of log(W / (H/3)) is shown in Fig. |221 According to the sim- 
ulation, the maximum- likelihood value of W{U(3) is 68 A, 
and the 90% likelihood interval for this quantity is as wide 
as 40-244 A. The observed value W(Hf3) = 330 ± 30 A is 
marginal, the likelihood of realizing a value greater than 
or equal to this one being «1%. This likelihood is of same 
order of magnitude as the 0.4% one to observe, as we do, 
no BSG and 2 WNL stars, given the mean IMF of the 
cluster. 

Using the same kind of simulations, we established 
that, for the same IMF slope and cutoffs, the fluctuation 
of W(H{3) would still be «20% for a cluster as massive as 
50 x 10 3 M Q . 

8. Conclusions 

In this work we have gathered the widest range of data 
available for the stellar cluster ionizing the giant HII re- 
gion NGC 588, both imaging and spectroscopy, covering a 
wide wavelength range from the ultraviolet to the far red. 

In the first part (Section [21 and 0}, we showed the im- 
portance of endeavoring to obtain and process as numer- 
ous and accurate observable data as possible, in order to 
constrain models of the most common form of OB asso- 
ciations within the limits of the relevance of these mod- 
els, rather than within the measurement errors and un- 
certainties. We discussed the problems encountered when 
fitting overall integrated nebular and stellar parameters, 
by means of a standard evolutionary population synthesis 
approach. In particular, by assuming an analytical IMF, 
the best fit model predicts the presence of BSGs simulta- 
neous to the required (since observed) existence of WNs, 
and consequently predicts too a low value of W(Hf3). This 
failure can be imputed mainly to two physical issues: the 
still important uncertainties in our knowled ge of the evo- 
lution of massive stars (as reviewed bv lMassevll2003|) . and 
the effects associate d with IMF sampling (as studied by 
ICervino et aLl f2002). We explored the latter in Section [5] 
through[7] In Scction|Sl we achieved quite a robust model 
of the cluster, and obtained good results. However, in this 
approach, we inferred a cluster age 50% higher than with 
the analytical approach, and a mass twice as large. In 
Section [7| we estimated the effects of departure of the 
observed set of stars from the one derived from the full 
sampling of the mean cluster IMF. In particular, we es- 
tablished that a diagnostic such as IF(H/3) is very sensitive 
to the BSG content of the cluster, and that the hardness of 
the ionizing radiation, an unobservable spectral property 
that plays an important role in photoionization models of 
nebulae, depends significantly on the WR content. More 
generally, we assessed that the radiation field of NGC 588 
is one among a large variety of possible SEDs for such a 
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cluster, and that this is easily explained by the following 
statements: the most massive stars dominate the flux of 
the cluster both in the observable and Lyman- continuum 
wavelength ranges, show a wide variety of spectral proper- 
ties, and are subject to strong population fluctuations. 

We suggest three practical considerations that will help 
to better understand and characterize the evolutionary 
stage of massive ionizing star clusters: (i) A "best fit 
model" statement may not be very meaningful per se, and 
should not make us complacent that it implies the best so- 
lution, (ii) When possible, use color-magnitude diagrams 
to obtain the integrated properties of the ionizing cluster, 
as opposed to assuming a perfectly sampled mass func- 
tion, (iii) Assess the importance of the few most massive 
stars on the overall SED of the cluster, both at ionizing 
and at non-ionizing wavelengths; this can be easily done 
by applying the lowest l umin osity limit criteria establishe d 
bv lCervifio et al.1 ll200.il) and lCervifio fc Luridianal l|2004h . 
If possible, use Monte-Carlo population models, provided 
they are relevant, and pending a more mature develop- 
ment of this approach to population synthesis. 

We also established that the initial mass distribution 
of the stars detected in NGC 588 is compatible with a 
Salpeter IMF, if it is treated as a stochastic process. 

Finally, in this work, we used the very common hypoth- 
esis of inst antaneous starburst. Howev er, recent works, like 
the one of lTenorio-Tagle et alJ l|2003j) . show that star for- 
mation in massive clusters spans time ranges such that 
the instantaneous formation hypothesis may not be valid 
for them. This statement does not call into question the 
significance of the disturbances caused by fluctuations in 
the high mass end of the IMF, but it is an additional un- 
certainty that is worth accounting for in further models of 
massive star clusters. 
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Appendix A: Magnitude interpolation 

At given age and metallicity, the isochrone data available 
from our programs consist in an array of physical param- 
eters of model stars, each line containing the following 
parameters for the current star: the initial mass M , the 
effective temperature T e g , the bolometric luminosity I/bob 
and the surface abundances, in particular the hydrogen 
one X. From these parameters, we wanted to derive ab- 
solute magnitudes, to be compared to the observations, 
from atmosphere models. In our case, f o ur kinds of mod- 
els we r e available: iLeieune et all l|l997|) . IPauldrach et alJ 
ll200lh . Iffillier fc Miller! lll998h and black body. For each 
of the first three kinds of models, whose selection for 
given physical parameters was operated the same way as 
Starburst99 does, stellar fluxes are tabulated for a fix ar- 
ray of wavelengths and for various effective temperatures; 
in the case of Lejeune and Pauldrach models, the sur- 
face gravity g is also present. When one of these grids is 
selected, Starburst99 performs nearest research interpo- 
lation along the logT off axis (for Hillicr grids) or in the 
(logT cf f, log 5) plane (for Lejeune and Pauldrach grids) to 
select a flux array, and normalizes the latter for its to- 
tal to be equal to Lboi- This non-continuous interpolation 
method can be exploited to compute the total spectrum 
of a model cluster, as the interpolation errors will tend to 
cancel each other. However, the resulting isochrone curve 
in a magnitude-magnitude of color-magnitude diagram 
will be irregular, with jumps resulting from the sudden 
transition from an array of the exploited grid to another. 
Here we suggest another form of interpolation to obtain 
more accurately determined model magnitudes. 

Let us first consider the case of the black-body emis- 
sion law: the wavelength-and-temperature-dependent flux 
is 

M\T) = 2 ^ eh J kT _ x (A.l) 



The change of variable x = hc/XkT = 1.44 x 10 8 /AT then 
gives: 

lnB A (A,T) = In (^^) ~ H e * ~ 1) (A.2) 

= ln(^-(x + ln(l-e-*)) (A.3) 

This means that at a given wavelength A, \nB\(X,T) is 
a linear function of the temperature-dependent variable 
y = \n{e x — 1) = x + ln(l — e~ x ). This is the main idea 
underlying the interpolation of magnitudes from the grids 
of model stellar spectra. Indeed, the wavelengths covered 
by the HST filters F547M, F439W, F336W and F170W 
are line-free and keep quite far from the continuum jumps 
(in particular, the Balmer jump); hence at these wave- 
lengths, we expect the stellar fluxes to behave nearly as in 
the black-body case. Consequently, we may calculate the 
flux predicted at wavelength A, for the effective tempera- 
ture T e g and the gravity g, with the fluxes available in the 
grid in use at the two temperatures T\ and Ti that "sur- 
round" T c ff, and at the same logarithmic surface gravity 
L = log g, using the following linear interpolation formula: 

mF A (A;T eff ,Z) = 
(y 2 - y) lnf A (A; T 1; L) + (y - m) lnf A (A; T 2 ,L) 
1)2 -yi 

This formula was straightly applied to the grids of Hillier 
model spectra, ignoring the L variable which is undefined 
for these grids. For Lejeune and Pauldrach grids, we also 
needed to interpolate the fluxes as functions of L. This 
was done the following way: 

— when, for the considered temperature, T\ or T2, at least 
two values of the surface gravity were available, the 
two ones closest to <?, namely g\ and 172, were used to 
perform linear interpolation as a function of L: 

In F x (A; T h L) = (A.5) 
(L 2 - L) In Fa (A; Tj,L{) + (L - L x ) In F A (A; Tj,L{) 
L2 — L\ 

— when only one value surface gravity was available in 
the grid (a marginal case that could occur with the 
Lejeune grids for the highest temperatures), its corre- 
sponding array of fluxes was adopted "as is" , i.e. ne- 
glecting the dependence of the stellar flux on gravity. 

When using Lejeune or Pauldrach grids, the interpolation 
algorithm was finally the following: 

1. in the grid, search for the two effective temperatures 
Ti and T 2 nearest the temperature of interest T e g ; 

2. for each of both T i; determine the two values gn and 
gi2 of the surface gravity closest to <?, and calculate 
the array of fluxes F\(\; Ti, log g) by linear interpo- 
lation in the "segment" [loggn, log ^2]- If actually 
only one value of the gravity gi is available, then set 
Fx{\]T i ,\ogg) = F x {X;T i ,\o E g i ); 
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3. at each wavelength, calculate y, y\ and 2/2, and apply 
formula (|A.4|I to determine the interpolated flux. Note 
that due to the fact that the variable x = 1.44 x 10 8 /AT 
could have values up to rj 800, we used the formula 
y = x + ln(l — e~ x ) and not y — ln(e x — 1), as the latter 
would have caused floating overflows when running our 
codes (e 800 - 10 347 ). 

As we saw above, our interpolation method is meaning- 
ful for wavelength ranges where the stellar fluxes behave 
nearly as black-body ones. However, it does not apply at 
wavelengths of lines or near discontinuities (e.g., Lyman 
and Balmer jumps). The interpolation proposed here is 
preferred to calculate smoother running isochrone magni- 
tudes, but does not involve a clear advantage to compute 
model spectra of clusters. 



